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ABSTRACT 

In this work, we study the formation and evolution of dark matter halos by means of 
the spherical infall model with shell-crossing. We present a framework to tackle this 
effect properly based on the numerical follow-up, with time, of that individual shell 
of matter that contains always the same fraction of mass with respect to the total 
mass. In this first step, we do not include angular momentum, velocity dispersion 
or triaxiality. Within this framework - named as the Spherical Shell Tracker (SST) 
- we investigate the dependence of the evolution of the halo with virial mass, with 
the adopted mass fraction of the shell, and for different cosmologies. We find that our 
results are very sensitive to a variation of the halo virial mass or the mass fraction of 
the shell that we consider. However, we obtain a negligible dependence on cosmology. 
Furthermore, we show that the effect of shell-crossing plays a crucial role in the way 
that the halo reaches the stabilization in radius and the virial equilibrium. We find 
that the values currently adopted in the literature for the actual density contrast at 
the moment of virialization, 5mr, may not be accurate enough. In this context, we 
stress the problems related to the definition of a virial mass and a virial radius for the 
halo. The question of whether the results found here may be obtained by tracking the 
shells with an analytic approximation remains to be explored. 

Key words: cosmology:theory — dark matter — large-scale structure of universe — 
methods : numerical 



1 INTRODUCTION 

In the hierarchical scenario for structure formation in the 
Universe, the small primordial density fluctuations grow due 
to non-linear gravitational evolution and finally become the 
first virialized structures (halos). In this picture, larger Cold 
Dark Matter (CDM) halos will be formed by the accretion 
and merger of those first smaller halos, forming in this way 
massive structures, and so on. This scenario, that constitutes 
the actual paradigm of hierarchical structure formation, is 
able to explain in general terms the universe that we see 
today. Yet, we do not have a framework or theory capable 
of reproducing this picture accurately. In this context, N- 
body cosmological simulations are a powerful tool to try 
to understand the formation and subsequent evolution of 
CDM halos. They constitute a very important help to build 
any theoretical model and their predictions explain many of 
different observations. 

Basically, there are two analytical approaches that make 
the problem tractable, although some simplifications have to 
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be made and, as it was said, comparison between these ana- 
lytical studies and simulations are crucial to make progress: 
the Press-Schechter formalism (Press & Schechter 1974), 
based on the role of mergers (Nusser & Sheth 1999; Man- 
rique et al. 2003), and the spherical infall model (SIM) fo- 
cused on the understanding of the collapse of individual 
objects. We must note that in the Press & Schechter for- 
malism, the SIM has also been widely used, but from a 
statistical point of view, to treat problems related to mass 
accretion histories, mass function, etc. The SIM, first devel- 
oped by Gunn & Gott (1972) and Gunn (1977), describes the 
collision-less collapse of a spherical perturbation in an ex- 
panding background. In those two articles, they introduced 
for the first time the cosmological expansion and the role of 
adiabatic invariance in the formation of individual objects. 
Fillmore & Goldreich (1984) and Bertschinger (1985) found 
analytical predictions for the density profiles of collapsed 
objects seeded by scale-free primordial perturbations in a 
fiat universe. Hoffman & Shaham (1985) generalized these 
solutions to realistic initial profiles in fiat and open Fried- 
mann models, and Baarden et al. (1986) (hereafter BBKS) 
improved this work introducing the peak formalism. Later, 
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some studies have been done to include more realistic dy- 
namics of the growth process of dark matter halos (e.g. Pad- 
manabhan 1996; Avila-Reese, Firmani & Hernandez 1998; 
Lokas 2000; Subramanian, Can & Ostriker 2000). 

In parallel, a large amount of numerical work have been 
done. Quinri, Salmon & Zurek (1986) and Frenk et al. (1988) 
obtained isothermal density profiles (p oc r~^) of CDM 
haloes, and Dubinski & Carlberg (1991) and Crone, Evrard 
& Richstone (1994) basically reproduced the predictions of 
Hoffman & Shaham (1985) and found some evidence for no 
pure power-law density profiles. Later, it was established 
that the density profiles of CDM halos have an universal 
form (Navarro, Prenk & White (1996, 1997, hereafter NEW), 
with p oc in the inner regions and p oc in the out- 
skirts, although there is still controversy about the shape of 
the profile near the center and recently it has been found 
that the profiles flatten out close to p oc beyond ~ 2 
virial radius. Moore et al. (1998,1999), amongst others, find 
p oc r^^'^ in the very center (although in a more recent 
work by Graham et al. (2006) they find p oc r^"'^ at 0.1 
kpc from the center) and other authors (Jing & Suto 2000; 
Klypin et al. 2001; Ricotti 2003) find an inner slope ranging 
from —1 to —1.5 depending on halo mass, merger history 
and substructure. Concerning the outskirts of dark matter 
halos, Prada et al. (2006) carried out a detailed study and 
concluded that a 3D Sersic three parameter approximation 
provides excellent density fits up to ~ 2 times the virial 
radius, although these profiles differ considerably from the 
NFW ones beyond 2 virial radius. 

There are also plenty of works in the literature using 
the SIM to predict the density profiles of dark matter halos 
mainly focused on explaining their central regions. More- 
over, the SIM has been widely used to obtain some quanti- 
ties specially relevant and directly related to crucial stages 
in the formation and evolution of CDM haloes for different 
cosmologies, redshifts, etc. A particularly relevant quantity 
is the value of the overdensity at the moment of virializa- 
tion Svir (A„jr usually in the literature), where overdensity 
is defined here as a number of times the background density, 
and its linear counterpart Si^^ir- The values of Si^^ir and Svir 
were obtained introducing the virial theorem into the SIM 
formalism. This has important implications in the way we 
define the virial radius (the radius that attains an overden- 
sity 5vir inside) of dark matter halos in N-body cosmological 
simulations. Svir is conventionally chosen to be near 180 for 
an Einstein-deSitter cosmology (e.g. Peebles 1980), or 340 
for the Q,A = 0.7 cosmology (e.g. Bryan & Norman 1998; 
Lokas & Hoffman 2000). 

In the standard derivation of Si^vir and Svir, the typical 
way to proceed is to assume that a shell of matter stabilizes 
at an epoch twice the time of turn-around (i.e. the time pre- 
dicted by the standard SIM to collapse into a point), and 
in average with a radius that is 1/2 the turn-around radius 
(e.g. Peebles 1980; Lacey & Cole 1994; Eke, Cole & Frenk 
1996). This 1/2 factor (in the Einstcin-deSitter cosmology; 
for other cosmologies we need to use the Lahav equation, 
Lahav et al. 1991) is called the collapse factor. However, the 
justification to introduce this collapse factor and to suppose 
the time of virialization as twice the time of turn-around, 
is poor and lack a solid theoretical background. In contrast, 
in this work we will study the spherical collapse without 
supposing any collapse factor, only taking into account the 



shell-crossing as the dominant effect. The angular momen- 
tum and velocity dispersion may also play an important role. 
The question is: if these effects were included in the model, 
would we obtain the same values for Svir and Si^vir that 
those found in the most simplistic scenario described by the 
standard SIM? This issue is one of the aims of the present 
work. 

The main goal of this line of work is to develop a theo- 
retical framework that help us to understand the dynamical 
elements that determine the process of formation of struc- 
tures (collapsed objects) using spherical symmetry to ex- 
plain main properties of dark matter halos. In this first work 
we will tackle these questions by means of a "cold" collapse, 
that is, without including the effects of the velocity disper- 
sion and angular momentum. The point is to ascertain if the 
non-uniformity of the density profiles generated via shell- 
crossing is able to provide the radial motions necessary to 
produce the virialization and stabilization in an appropri- 
ate time scale. In a future work, we will include the angular 
momentum and velocity dispersion to go a step further. 

There are some issues that it is worth mentioning and 
that make this work different from previous works that also 
included the shell- crossing in their formalism (e.g. Lokas 
& Hoffman 2000; Nusser 2001; HioteUs 2002; Ascasibar et 
al. 2004) . The way to proceed in these works is to handle the 
effect of shell-crossing by means of an adiabatic invariant, 
once the standard SIM becomes incorrect for late stages of 
the evolution. This adiabatic invariant, also known as ra- 
dial action, makes the problem analytically tractable, and is 
based on the fact that the potential evolves in a time larger 
than the orbital period of the most inner shells. In contrast, 
we will study the shell-crossing effect doing a follow-up of 
the radius that contains inside always the same fraction of 
the virial mass. This way to tackle the problem is only one 
of the possible options, but is essential, for example, in order 
to build and study the relationship between the actual en- 
closed density contrast S, defined as S — p(<'^^~<^"'> ^ with 
< Pm > the mean matter density of the Universe, and the 
linear density contrast. Si, obtained from the linear theory. 
Only Gehard Lemson did something similar, although using 
N-body simulations and mainly focused on showing how ac- 
curate arc the predictions of the standard SIM compared 
to his simulations (Lemson 1995). Despite the fact that he 
showed that the SIM is a powerful tool to understand the 
evolution of halos, he never provided detailed quantities and 
relations for the actual and linear density contrasts. The 
function 5i{S) is very important to obtain the density pro- 
files of dark matter halos, as we discussed in previous works 
(Prada et al. 2006; Betancort-Rijo et al. 2006). Sheth & Tor- 
men (2002) parametrized this function for the standard SIM. 
The framework presented here will allow us in the near fu- 
ture to provide also a simple parametrization for Si{S), but 
taking into account the important effect of shell-crossing, 
together with others relevant effects such as the angular 
momentum and velocity dispersion. This will lead us, for 
example, to obtain Si^vir and its corresponding S^ir, to ex- 
plain the shape of the dark matter density profiles or to 
shed light on the mass functions. All of this without sup- 
posing any collapse factor, as pointed before, or other vague 
assumptions. Nevertheless, it will not be possible to obtain 
useful applications for the moment, since in this first work 
we will include in our study only the shell-crossing, which 
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is the dominant effect. The full treatment will be done and 
presented in an upcoming work. Here we will provide the 
first results of our theoretical framework related to the role 
played by the shell-crossing. 

The paper is organized as follows. In section[2]we briefly 
describe the SIM and explain the formalism and unities that 
we will use in the rest of the work. In section [3] we will study 
in detail the dependence of the way that the evolution oc- 
curs varying some parameters, in particular, the virial mass 
of the halo, the fraction of mass for a given virial mass, 
and the cosmology. Section |4] will be specially dedicated to 
the moments of virialization and stabilization according to 
a given criterion, and their dependence with the same pa- 
rameters described above. We will also emphasize here the 
difference between these two concepts. Finally, in section O 
we address the main results and ideas of the work, and point 
the lines for a future work. 



2 THE SPHERICAL SHELL TRACKER 
FRAMEWORK 

In this section, we first describe the standard SIM and its 
equations, and then we present the formalism that we will 
use in the rest of the work, which will imply to describe 
the density profile and define our own units. This will allow 
us to handle easily the equations involved. Later, the algo- 
rithm that we used to obtain the results will be described 
carefully step by step. All together will be known as the 
Spherical Shell Tracker Framework (SST). The main objec- 
tive of this section is to make easier a possible reproduction 
and implementation of the SST framework. 



2.1 The formalism 

In a flat Universe with SIa = 0, the evolution of a homoge- 
neous spherical (positive) density perturbation (the simplest 
way to tackle the problem of structure formation) with a 
mass M and radius R, is given by Newtonian dynamics (as 
shown by Tolman 1934 and Bondi 1947), provided that R 
be much smaller than the Hubble radius: 



<fR 
tain; 



-G M 



(1) 



Integrating, since M is constant by definition, we ob- 



1 (dRV 

2 \ dt) 



G M 
R 



= E 



(2) 



where E determines whether the sphere expands forever 
{E > 0) or it finally contracts {E <Q). 

We can describe in more detail this spherical perturba- 
tion with a large number of mass particles, and even it is 
possible and more useful to imagine these particles as con- 
centric shells (thanks to spherical symmetry) that do not 
cross each other, and each of them with a radius r{j,t), 
where j denotes the shell, which satisfies equation ((T|: 



d^r{j,t) _ -G M{j,t) 



(3) 



where M{j, t) is the enclosed mass for each shell j at 
time t: 



M{j, t) = p^ru (|^r(j, t)3) [1 + S{r{j, t))] (4) 

being pcrit the critical density of the Universe, and S{r) 
the actual density contrast within r{j,t): 



pcrit 



3 H\ 
Stt G' 



S{r) 



p{< r)- < Pm > 

< Pni > 



(5) 



with H the Hubble constant, and < p,n > the mean 
matter density of the Universe. 

As long as shell- crossing does not occurs, the actual 
density contrast is related to the linear one (given by the 
linear theory, see e.g. Padmanabhan 1993) in the Einstein- 
deSitter cosmology by the formula (Sheth & Tormen 2002): 



1.68647- 



1.35 



1.12431 



0.78785 



(6) 



The inverse function, 
al. 2004); 

5{Si) = 0.993[(1 - 0.607(5; - 
+ e{S, - 1.55))Sf)r^-'''' ~ 1 
being 9 the step function; 
1 if X > 



(l+5)2/« (l + <5)l/2 (1+5)0.58661 

S(5i), is given by (Patiri et 



6.5 X 10^^(1 -6'((5;) + 



(7) 



0(x) = 



if X < 



It is possible to make some simplifications in the equa- 
tions, choosing in an appropriate manner the value of some 
parameters. In particular, we choose; 

time unit = initial time 
length unit = initial radius of the protohalo, Ri 
mass unit = [1 -I- S{Ri)] 

According to these units, and taking into account equa- 
tions @ and ([SJ, and an Einstein-deSitter cosmology (where 



H 



It- 



we have; 



Pcrit. i — . , 

47r 



r- 2 



(8) 



where i refers to the initial time. 

The Lagrangian radius q for each shell j (i.e. the co- 
moving radius at t ^ 0) is related to the Eulerian one r 
by; 



q{j)=nij) [l + 5(r,(j))]- 



(9) 



So, for the initial enclosed mass of a shell j, we now 
have simply (in our units and using Eq.(|4)l); 

M{j,U)=M{j)=qijf (10) 

We must note that this enclosed mass of a shell j, M{j), 
is different from the mass of each shell, Msheu{j)'- 



MsHellij)=M{j)-M{j-l) 

for r{i) ^ r{j) 



M{j,t) ^y^Msheiiii) 



1=1 



(11) 



whit n the total number of shells. 

To obtain q{j) using Eq.([9]) we need S^{ri{j)), that is, 
the actual density contrast at initial time, and to this end we 
need the linear profile at initial time, Sl(q{j)). In this work 
we will use the linear profile presented in Prada et al. (2006) 
and Betancort-Rijo et al. (2006); 
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Table 1. Values of b and Q necessary to use the approximation 
for 5i{q) given by l|14| l. 



M (/i-1Mq) 


Q (h-iMpc) 


b 


6.5 X IQi" 


0.57 


0.1889 


5 X 10" 


1.1252 


0.2202 


3 X 10^2 


2.0445 


0.2544 


2 X 1013 


3.848 


0.301 


5 X 10^"' 


11.125 


0.41 



where Zi is the redshift at initial time. We can obtain 
5^{ri{j)) from 5l{q{j)) using the function given in ([7]). In- 
serting this S{5l{q{j))) in Eq.® we obtain the Lagrangian 
radius for each shell j, q{j), and also using Ea. (|10|) its en- 
closed mass M{j). 

Once we have the expressions related to the initial 
conditions and we have presented the density profile, we 
need the equations of evolution. If the shells do not cross 
each other, then there is an analytical solution for ((2)1 (e.g. 
Martinez & Saar 2002) that can be written in the parametric 
form: 



5l{q) = Sl^yir 



"•12(g) 



(12) 



r = rc(l — coari); t — tdj] — sinri) 



where 5i,vir is the linear density contrast at the moment 
of virialization, q and Q are the Lagrangian radii related to 
r and Rvir respectively, that can be obtained using equation 
|[9}, and: 



{cj{x)f = ^ / 1 I' w^[xk) e dk 



<Ti2 = <Ti2(g) = 2^ I 1 Wriqk) WriQk) k^ dk 

, , 3(sm X — X cos X ) 
Wt(x) = ^ 3 - 



(13) 



where stands for the power spectra of the density 
fluctuations linearly extrapolated to the present. 
There is a good approximation for 5i{q): 



Si{q) = Si^yir exp 



_1 dlna(x) 
2 dlnx 



(14) 



with b a constant depending on the mass. In Table \T\ 
we present the values of b and Q that we use for each mass. 
Moreover, it is necessary to assign a value to Si^^ir so that we 
can use the density profile, although one of our final goals 
is to obtain a precise value for it. In this work we used a 
Si.vir ~ 1.9, a value which led to good results in previous 
works (Prada et al. 2006; Betancort-Rijo et al. 2006). 

Essentially, the profile given in Ea. (ll2|l and its approx- 
imation in Eq. (|14|l takes into account only the restriction 
SiiQ) = Si^^ir. In Hiotelis (2002) and Ascasibar et al. (2004), 
a very similar density profile was also used, but using the 
BBKS peak formalism to compute the initial conditions. In 
a future work we will use a more sophisticated density pro- 
file that includes also the restriction Si{q) < Si^^ir for q> Q 
(see Betancort-Rijo et al. (2006) for a more detailed descrip- 
tion), resulting in steeper actual density profiles for smaller 
masses (as confirmed by numerical simulations, see Prada 
et al. 2006). This fact will probably change slightly the re- 
sults. For the sake of simplicity we preferred to use a simple 
profile now, although the major results of this work will not 
depend on the assumed profile. 

For the Einstein-deSitter cosmology, the 5; (g(j)) profile 
can be obtained from equation p2p simply rescaling by: 



SKqm = 



1 



1 + 



(15) 



where: 



GM 



f-c — 1 



dt 
drj 



R 

c 



(16) 



(17) 



Here c is the velocity of hght and there is a change to 
a non-dimensional variable r?. This solution means that the 
shell expands until it reaches a maximum radius rta, the 
turn-around radius, at a given time tta, which is different 
for each shell, and after that point the shell starts to con- 
tract. We can integrate analytically equation ([3]) to study 
the evolution of the spherical density perturbation, at least 
until the turn-around, thanks to the fact that the enclosed 
mass of the shells do not change with time. Nevertheless, 
after the turn-around, the re-collapse begins and the shell- 
crossing also starts, so we can not proceed in the same way. 
At that point, it is common to use a prescription based on 
an adiabatic invariant, to account for this secondary infall 
and shell-crossing (e.g. Lokas & Hoffman 2000; Nusser 2001; 
Hiotelis 2002; Ascasibar et al. 2004). On the other hand, one 
can also integrate numerically the equation computing 
at each time step the new radius, velocity and enclosed mass 
for each shell. This is the method that we use in our work. 
Our purpose is to study and to include the shell-crossing in 
our treatment in a natural way, i.e. without making any as- 
sumption about the collapse factor, the time at which virial- 
ization occurs, or any other simplification or approximation. 

We first divide our spherical density perturbation in n 
equal spherical shells (equivalent to particles), all of them 
with the same thickness, and later we choose the shell j that 
contain a given fraction of mass of the total protohalo. We 
do so for every time step, from the start of the evolution 
to the end: we recompute the new enclosed mass for each 
shell at each time step, and we always select that one that 
contains the fraction of mass we are interested in (in that 
sense, n must be big enough to choose with high precision 
and without problems at each step a shell that contains ex- 
actly the required fraction of mass; in our case, n = 3000 
was enough). If we follow for a long time the shell related to 
this fraction of mass, at the end its radius will be almost con- 
stant (although the corresponding physical shell will change 
with time), that is, we will reach stabilization (see section 
4). Lu et al. (2006) used a similar algorithm, but they di- 
vided the halo in equal mass shells, instead of shells with the 
same thickness, as we do. Moreover, their motivations were 
different, mainly focused on explain the inner shape of the 
density profiles, and they did not carried out a follow-up of 
any shell in particular. 
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2.2 The algorithm 

We now describe the algorithm we used to compute the rele- 
vant quantities (radius, velocity and enclosed mass) for each 
shell at each time step. 

First, we need to obtain the initial conditions: 

(i) We divide the protohalo in n equal shells to calculate 
our array of initial radii, ri{j), that contains the radius r 
for all the shells. Remember that, in our units, the radius 
of the total cloud is unity; moreover, j increases decreasing 
the radius, so r{j = 1) = Ri (the radius of the whole halo), 
and r{j — 3000) is the radius of the deepest shell. 

(ii) In a first approximation, we make q{j) = ri{j). This 
will allow us to compute a second and better estimation 
for q(j) using in an appropriate way the relation given by 
equation ((Ojl, that is: 



(18) 



where we introduce in the right side the q{j) as given by 
the first approximation. The function 5{5i) is given by Eq.([7]) 
and 5l{q{j)) is given by (fT5|) . 

(iii) Now, we will use the q{j) obtained in the last step as 
a new approximation to compute again a better estimation 
for q{j), according to ea. (|18p . 

(iv) Step (iii) must be repeated until there is no difference 
between the q{j) that we obtain after each iteration, or this 
difference is less than at least 5% between two consecutive 
iterations. 

(v) With ri{j) and the last and best estimations for the 
initial Lagrangian radii of the shells, q{j), we can calcu- 
late the initial enclosed mass array, M{j), using equation 
HlOp. and the mass of each shell, Msheii{j), knowing that 
MshMj) = M{j) - M{j - 1). We will need MshsiiU) later 
to compute the enclosed mass array at each time step, since 
the mass of each shell will be always the same, although the 
order of the shells will be modified. 

(vi) We also need the initial velocity for each shell, Vi{j), 
as given by the SIM (Betancort-Rijo et al. 2006): 



31+S{5l{q{j))) 



(19) 



where 5i{5) and S{5i) are given by Eq.® and Eq.Q, and 
Sl{q{j)) is given by Ea. (|15|l . Hi is the Hubble constant at 
initial time, which in our units is Hi = 2/3. 

(vii) Now it is possible to select the shell that contains 
the fraction of mass that we are interested in. Our studies 
will be focused on the follow-up of this shell in particular. 

Once we have calculated the initial conditions, we now need 
to obtain the equations of evolution: 

(viii) In our units, equation (|3} can be written as: 



dV(j,t) _ 2M{j,t) 



(20) 



dt^ 9r{j,t)^ 
so the equations of the evolution for r{j) and v{j) are: 
r{j, t + At)= r{j, t) + v{j, t) At (21) 



/ • X 2M(j,t) At 
v{j,t + At)^v{j,t)~-- ' 



(22) 



9 r(j,t)2 

(ix) We also need to compute how the linear and actual 
density contrasts evolve with time: 



Sliht) - SliqUY) t^ 



[l + 5{Sliqm] 



(23) 
(24) 



(x) At each time step, we need to recalculate the new 
enclosed mass array, M{j), using the shell mass array 
Msheiiij), since the enclosed mass for a given shell j is equal 
to: 



M(j,i) 



sheinl) 



for r{i) !^ r{j) 



(25) 



(xi) At this point, we can select again the shell that con- 
tains that fraction of mass we want to study, to see what 
happens with its radius, velocity, and linear and actual den- 
sity contrasts. 

(xii) For each time step, we will have to repeat (ix) to 
(xii). 

It is worth mentioning, for possible reproductions of the 
results, that we used an optimized temporal step of 0.003 
in our units, which is good enough to give us robust results 
of Si and 5 (we checked these values using well known mo- 
ments of the evolution like the turn around, where there is 
no shell-crossing yet). Moreover, the beginning was set to an 
initial redshift Zi = 15, to make sure that we are still well 
inside the linear regime, i.e. the initial value of 5i is small 
enough. However, it must be noted that by z = we do not 
necessarily mean the present time. When we are considering 
a given mass scale, z = corresponds to the time of viri- 
alization of that scale, that is, the time when Si within the 
Lagrangian virial radius is ~ 1.7 (more precisely, our z — 
is that one where we have a linear density profile given by 
Eq. ll4|) . So, at z=l/15. Si in that scale is much smaller than 
one. 

Furthermore, there is another important question that it 
is necessary to take into account to implement without prob- 
lems the described algorithm. This is the fact that we have 
not included yet the effect of the velocity dispersion and 
angular momentum. Therefore, we are in a totally radial 
(cold) collapse and we will have problems in the very center 
of the halo if we simply integrate numerically the equations 
following this framework. When we compute, according to 
equations (|2ip and (|22[). the new radius and velocity of a 
shell which is located very near to the center, we can obtain 
at the following time step a negative radius and a positive 
velocity, which means that actually this shell has crossed 
through the center and now it goes from the inner regions 
of the halo to the outer ones. In these circumstances, en- 
ergy is not exactly conserved due to numerical reasons. The 
distance that the shell covers in only one time step is com- 
parable to its radius, which gives a considerable "leak" of 
energy. There are different ways to solve this problem; one 
of them, the solution we chose, is to define a parameter m 
to measure properly this effect and to help us to prevent 
this lost of energy. If we define for each shell the parameter 
m as m = Ar/r{j,t) (where Ar = r{j,t + At) — r{j,t) is 
the distance that the shell has covered along this time step). 
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then there is no problem while m is small enough, but when 
m reaches a larger value, the way to minimize the lost of 
energy in the process is to change the velocity by its ab- 
solute value, and keeping intact the value that we have for 
the radius. Doing so, we skip the very center and the loss of 
energy will be minimum. After many attempts, we saw that 
a value of m = 0.02 leads to very good results. Once ve- 
locity dispersion and angular momentum are included, this 
parameter m will not be necessary. 

In the AC DM cosmology, the formalism and the algo- 
rithm are the same, but we must introduce some changes in 
the initial conditions and in the equations of the evolution 
the spherical to take into account the different cosmology 
with A 7^ 0. The equations, modified adequately, are pre- 
sented in Appendix A. 




0.2 



0^~w \ ' >: n 'f 



3 THE EVOLUTION OF THE HALO: EFFECT 
OF SHELL-CROSSING 

At the beginning of the evolution, we will not obtain any 
difference using our formalism or using the standard SIM, 
because there is no shell-crossing yet. However, when this 
effect starts it is clear that this will become false. However, 
the expected deviation with respect to that given by the 
standard SIM may be different if one uses different values 
for the virial mass of the halo, or study different fractions 
of mass with respect to this virial mass, or if we move to a 
different cosmology. To this end, that is, to quantify in detail 
how big are the dependences on these factors, we carried out 
a study varying the virial mass of the halo, M^ir, the fraction 
of mass related to this virial mass, that we call Mjrac, and 
the cosmology through the value of Q.a- 

To illustrate the way in which shell-crossing occurs and 
affects to the evolution of the halo, in Figure [1] we show the 
evolution with time of the radius related to different Mfrac 
for the particular case of an Einstein-deSitter cosmology and 
a virial mass of M^ir = 3 x lQ^^h~^ Mq. Both the radius and 
time are expressed in units of the turnaround values (so we 
can compare between different Mfrac in the same scale). It is 
worth mentioning that Lemson (1995) presented in the same 
way data from his simulations and he obtained essentially 
the same results as shown here for the evolution of individual 
shells. 

Before the first shell-crossing happens, the behaviour 
of the radii of different Mfrac is essentially the same. This 
first shell-crossing occurs just before twice the time of turn- 
around (the time of virialization for the usual models), and 
what we can see is that this first shell-crossing means the 
beginning of the stabilization in radius, which finally occurs 
some time after that (in next section, we will carry out a de- 
tailed study of this process together with the virialization). 
The larger radius oscillations for each curve beyond ~ 'it/ta 
(see Figure 1) are only noise due to the growth of numer- 
ical errors with time, although in the case of Mfrac = 1 
the larger deviations at larger times are partially and prob- 
ably due to border effects, i.e. the shell that contains the 
required fraction of mass is near the border of the halo at 
that time, so there are no enough shells above to obtain a 
good behaviour using our algorithm. 

In Table [2] we summarize the results for the linear and 
actual density contrasts for a halo with virial mass M^ir = 



1 2 3 '1 5 

t / t,. 

Figure 1. Evolution with time of the radius for different Mfracy 
for a halo with a virial mass M^ir = 3 X IO^^/i^^Mq and an 
Einstein-deSitter cosmology. Both radius and time are in units of 
the turnaround radius and time respectively. From down to top, 
the curves are for Mfrac = 0-2; 0.5; 0.8; 1.0. Triangle means 
i5 = 180, squares the first shell-crossing, and circles indicate the 
time of collapse according to the standard SIM. The horizontal 
dashed-line corresponds to half the turnaround radius, and the 
vertical one the time of collapse, i.e. twice the turnaround time. 

3 X W^^ Mq and five values of Mfrac for two different 
cosmologies: the Eintein-deSitter and a model with SIa 7^ 
0. In each case, the corresponding values of Si and 5 are 
given for critical or interesting moments of the evolution, in 
particular when the first shell-crossing occurs (FSC), when 
S = 180 (A180), S = 340 (A340) and when the collapse 
occurs {2TA) according to the standard SIM, that is, twice 
the time of turn-around. The selection of A180 and A340 
was done because they are the preferred values of S^ir in 
the literature for a flat universe with 57a = and Qa ~ 0.7 
respectively. They would be also useful to show a possible 
dependence (or not) of the function 5{Si) on the cosmology. 
We must note that, although Table[2]is only for a given virial 
mass, the same kind of study was done for the evolution 
of halos with virial masses M^ir = 6.5 X 10^°/i"^Mq and 
Mvir — 5 X 10^''/i~^Mq. We observed the same tendencies 
in the data and achieved the same conclusions. In Table [3] 
we show the values of Si and S that we obtained for the 
three mentioned virial masses and for the Eintein-deSitter 
and the SIa 7^ cosmologies, for the particular case in which 
we fixed AIfrac=0.5. 

It is worth mentioning that we did all the calculations 
to get a cosmology with Qa = 0.7 at our nominal z = 0, so 
in general the values of Si and 5 given in the tables are actu- 
ally related to other values of the cosmological constant, i.e. 
the value of this constant at the corresponding time FSC, 
A180, A340 or 2TA. Because of this fact we give in the 
tables the corresponding value of /3 = Qm/i^A at the cor- 
responding moment of evolution in the case of a Qa 7^ 
cosmology (in the Einstein-deSitter case it is not necessary) . 
From the value of /3 we can easily deduce also the redshift 
at which that moment of evolution occurs, knowing that 
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P = Vlm/^A = Po (1 + z)^, with /?o = f^m.o/f^A.o IS the 
value of P at the time corresponding to our z = (which is 
Po = 0.429 for the n,„ = 0.3, = 0.7 cosmology). Higher 
values for P compared to this Po correspond to a virializa- 
tion or stabilization that occurs at earlier times {z > 0), 
while lower values represent future times (2 < 0). 

Some interesting conclusions can be inferred from the 
Figure 1 and Tables [2] and |3] The most important ones can 
be summarized as follows: 

(i) FSC , the first shell-crossing, occurs earlier as Mjrac 
increases. 

(ii) FSC occurs also earlier for larger masses. 

(iii) FSC sets the beginning of the stabilization in radius. 

(iv) FSC always occurs after A180 and A340 but always 
before 2TA, independently of M^ir and Mjrac- 

(v) A180 and A340 have essentially the same associated 
linear density contrasts, no matter the value of Mjrac or 
Mmr- This is because in all the cases, there has not been 
any shell-crossing before A180 and A340, so the standard 
SIM is still vaUd. 

(vi) Concerning the linear and actual density contrasts for 
a given Mfrac and M^ir, there is no substantial difference 
between the values obtained for different cosmologies. 

Specially relevant is the last conclusion, which means that 
there is no dependence with the cosmology in the values of 
5i and 5, or this dependence is really small and negligible. 

Furthermore, we must note here the very high values 
found for the actual density contrast 5 at the moment of the 
first shell-crossing {FSC) and collapse {2TA). The reason 
for that is that there are other important effects, together 
with shell-crossing, involved in the formation and evolution 
of dark matter halos and that we have not included yet in 
our model. In particular, angular momentum and velocity 
dispersion will become very relevant and by sure will reduce 
the values that we obtain for S. In fact, Avila- Reese, Firmani 
& Hernandez (1998), Hiotelis (2002), Ascasibar et al. (2004) 
and Shapiro et al. (2004), amongst others, introduce and 
study the angular momentum and find shallower density 
profiles in the inner regions, as expected. Hence, it will be 
absolutely necessary to take into account at least these two 
effects if we want to go a step further in our analysis and 
if we want to obtain a good and accurate parametrization 
for the function &i{5). Nevertheless, the framework and algo- 
rithm we are using, as well as the conclusions and tendencies 
we can obtain only including the shell-crossing, are totally 
valid although we can not reach, by now, exact values. In- 
cluding other effects in our framework, specially those ones 
mentioned above, will be part of a future work. 

Then, for the moment, we will not be able to pro- 
vide an exact relation between the linear and actual den- 
sity contrasts, neither a parametric form for the function 
Si(5). However, we can have a look to the relation that we 
obtain at this moment between both density contrasts, and 
try to extract some conclusions. In Figure [2l the function 
Si{5) is represented for the three virial masses under study 
and for the Einstein-deSitter cosmology. Figure[3]represents 
the same function but for different cosmologies, in partic- 
ular for the Einstein-deSitter case and SIa 7^ 0. The linear 
region is clearly visible in both figures below 5; ~ 1. In this 
regime there is no still any difference between the different 
curves and the value of 5 grows very slowly with 5i, as ex- 




0,0 0.5 1.0 1.5 2.Q 2.5 



Figure 2. The relation (5; —5 for three virial masses. From down to 
top, the curves correspond to M^ir = 5 X lCi^'^h~^MQ, M^ir = 
3 X IQi^h-iMQ and M„i^ = 6.5 X lO^'/i^iM© (an Einstein- 
deSitter universe and M frac=0.5 was used in all the cases). 




O.Q 0.5 1.0 1,5 2.0 2.5 



Figure 3. The relation 5i —S for two different cosmologies. From 
down to top, the curves correspond to the Einstein-deSitter case 
and to the 7^ cosmology (a virial mass of M^ir = 3 X 
W^^h~^ Mq and M frac=0-5 was used in all the cases). 



pected. Then, there is a phase where 6 increases very fast 
for small differences in 5i, starting from Si ~ 1.6 in all the 
cases. From this moment, the dependence with virial mass 
becomes clearly visible in figure (2] where we observe that 
the smaller the mass, the larger the values of 5 attain for 
the same value of 5i . Respect to the dependence on different 
cosmologies, it seems clear (see Figure |3]) that this depen- 
dence is really small, as already mentioned. 
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Table 2. Linear (Si) and actual (5) density contrast values related to some important moments in the evolution of a halo with a virial 
mass Myir = 3 X 10^^h~^ Mq and for the Einstein-deSitter and £7^ ^ cosmologies. /3 stands for the ratio between Qrn and Q/^ at the 
corresponding moment. See text for details, page 6. 



Einstein-deSitter (Qm = 1, f^A = 0) 



Moment 


Mfrac 

Si 


= 0.2 

5 


Mfrac 

Si 


= 0.5 

5 


-^^ frac 

Si 


= 0.8 

5 


Si 


= 1.0 

5 


Mfrac 

Si 


= 1.3 

5 


FSC 


1.667 


1840 


1.652 


770 


1.641 


500 


1.634 


398 


1.628 


329 


A180 


1.601 


180 


1.602 


180 


1.602 


180 


1.602 


180 


1.603 


180 


A340 


1.628 


340 


1.628 


340 


1.629 


340 


1.628 


340 


1.628 


340 


2TA 


1.695 


1426 


1.695 


653 


1.696 


424 


1.696 


350 


1.696 


289 





Mf 


rac 


0.2 


Mf. 


rac 


0.5 


Mf. 


rac 


0.8 


Mf 


7'ac 


1.0 


Mf 


rac 


1.3 


Moment 


Si 


5 


/3 


Si 


5 


/3 


Si 


5 


/3 


Si 


5 


/3 


Si 


s 


/3 


FSC 


1.670 


1728 


1.686 


1.655 


761 


1.305 


1.644 


516 


1.048 


1.639 


427 


0.910 


1.633 


356 


0.739 


A180 


1.606 


180 


1.958 


1.606 


180 


1.477 


1.605 


180 


1.162 


1.605 


180 


1.000 


1.605 


180 


0.800 


A340 


1.632 


340 


1.844 


1.632 


340 


1.382 


1.631 


340 


1.084 


1.631 


340 


0.929 


1.631 


340 


0.741 


2TA 


1.698 


1364 


1.582 


1.698 


632 


1.175 


1.697 


446 


0.912 


1.696 


380 


0.777 


1.696 


320 


0.612 



4 STABILIZATION AND VIRIALIZATION 

In the standard SIM and an Einstein-deSitter cosmology, 
the value of 5i corresponding to the final stage of evo- 
lution, to the so-called vtrialization, is usually taken as 
Si,vir ~ 1.686, that corresponds to an actual density contrast 
Syir ~ 180 (e.g. Peebles 1980). For the D,m = 0.3, ^a = 0.7 
cosmology, Si^^ir — 1.676 and S^ir ~ 340 (e.g. Lacey & Cole 
1994; Eke, Cole & Frenk 1996; Bryan & Norman 1998). As 
pointed in previous sections, those calculations are based 
mainly on the following assumptions: 

(i) The halo virializes within a radius that is, on average, 
a given fraction of its maximum radius (the turnaround ra- 
dius). This fraction is the collapse factor, and is equal to 
1/2 in the Einstein-deSitter cosmology, and in other cos- 
mologies it can be inferred from the Lahav equation (Lahav 
et al. 1991). 

(ii) The time at which virialization occurs is twice the 
time of turn-around, that is, the time at which the collapse 
happens according to the standard SIM. 

Although the values inferred for Si^rir and S^jir in this 
way are commonly accepted as the correct ones and are 
widely used in the entire literature, the reasons to make 
the assumptions given above lack a solid theoretical base 
(see section O. In fact, there are some works that point 
to another direction and estimate other values of 5i,„ir and 
S^jir- Jenkins et al. (2001), for example, find a better agree- 
ment with the simulations if S^ir is taken constant for all 
the cosmologies and near the value that it takes in the Ein- 
stein deSitter cosmology [S^ir ~ 180). Also Avila- Reese, Fir- 
mani & Hernandez (1998) find that a different value of Si^^ir 
with respect to those obtained using the above assumptions 
makes better the comparison between the analytical Press- 
Schechter mass distribution and the results of N-body sim- 
ulations. 

Moreover, there is another important question related 
to the virialization that should be considered here. In the 



framework of the standard SIM, it is possible to apply the 
virial theorem if we suppose the halo to be an isolated sys- 
tem. However, real halos are non-isolated systems, with sur- 
rounding material continuously falling or escaping from the 
system. Hence, the virial theorem at least in the standard 
form could not be applied in this case. Despite of this fact, 
the standard SIM together with the virial theorem have 
been used to obtain the values of Si^vir and 5vir, and these 
values have been taken as the references to define the virial 
radius and the virial mass of the halos, which is specially 
adopted in N-body simulations. Furthermore, this fact has 
been traditionally supported for radial velocity early studies 
of massive dark matter halos from simulations (Crone et al. 
1994; Cole & Lacey 1996) . These studies apparently showed 
that the virial radius in this way defined (using Si^^ir = 1.69 
and Sjjir ~ 180 in the Einstein-deSitter) constitutes an ad- 
equate boundary to separate the inner region of the halo 
in dynamical equilibrium, i.e. that region where the radial 
velocities are zero, from the external region showing infall 
velocities. The popularization of these ideas came contem- 
poraneously with works that defined the virial mass and 
virial radius in simulations according to these preliminary 
results (specially since the NFW papers). But the fact is 
that, as recently shown in Prada et al. (2006), this may not 
be totally correct. Concerning galaxy-size halos, for exam- 
ple, they display all the properties of relaxed objects up to 
~ 3 virial radius and there is no indication of infall of ma- 
terial beyond. Therefore, there is no reason to believe that 
only inside the virial radius, as currently defined, the halo is 
in equilibrium. In this context, it is important to understand 
the process of virialization more in depth. 

In our work, no assumption is done related to the virial- 
ization. No collapse factor, no time for virialization a priori 
is adopted. Including shell-crossing in the way we do will al- 
low us to obtain Si^vir and its corresponding Svir in a natural 
way, i.e. studying the evolution of different shells of the halo 
according to the SST framework, presented in section 2. It 
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Table 3. Linear (5;) and actual (5) density contrast values related to some important moments in the evolution of a halo for three 
different virial masses (M^i,- = 6.5 X lO^''?t~-'^M0, Alyir = 3 X 10^^ Mq and Myir = 5 X 10^'^h~^ AIq) and two different cosmologies 
(Einstein-deSitter and S7a 0). A value of Mf^-^c = 0.5 was set in all the cases. /3 stands for the ratio between Cm and Ha at the 
corresponding moment. See text for details, page 6. 

Einstein-deSitter (Hm = 1, f^A = 0) 
M„ir = 6.5 X IQiOft-^Afo A/i,ir = 3 X lO^^ft-U/Q Af„i^ = 5 x lO^/i-^Afo 



Moment 


&i 


<5 


ii 


<5 




<5 


FSC 


1.660 


1203 


1.652 


770 


1.631 


377 


A180 


1.602 


180 


1.602 


180 


1.601 


180 


A340 


1.628 


340 


1.629 


340 


1.628 


340 


2TA 


1.695 


986 


1.695 


653 


1.694 


335 



Ca 7^ 







6.5 X 10^° 


H-^Mq 




3 X 10^2 






5 X 10" 


ft-lAfg 


Moment 




<5 


/3 




<5 


/3 


Si 


.5 


/9 


FSC 


1.663 


1170 


1.160 


1.655 


761 


1.305 


1.637 


392 


1.718 


A180 


1.606 


180 


1.341 


1.606 


180 


1.477 


1.606 


180 


1.848 


A340 


1.632 


340 


1.254 


1.632 


340 


1.382 


1.632 


340 


1.737 


2TA 


1.697 


978 


1.063 


1.698 


632 


1.175 


1.698 


338 


1.489 



must be noted, however, that these values are still prelimi- 
nary, since it will be necessary to include the effects of angu- 
lar momentum, velocity dispersion and triaxiality to obtain 
precise and useful values. Nevertheless, this study will be 
suitable to isolate the role of shell-crossing and will be able 
to extract important conclusions related to the stabilization 
and virialization. We believe that these conclusions will not 
change when we introduce other physical considerations into 
the framework. 

It is important to note here the difference between these 
two concepts: stabilization and virialization. The first one 
can be inferred studying the behaviour of the radius of a 
given shell that contains a given fraction of the virial mass 
with time, as was shown in Figure [1] A criterion must be 
imposed to say if a given shell reaches stabilization or not, 
and when. The second concept, the virialization, will have 
to be inferred according to the virial theorem. There is no 
reason why stabilization and virialization should coincide, 
although instinctively one expect that they should be near 
in time at least. 



4.1 Stabilization 

In first place, we define the time of stabilization as the time 
at which the radius of the shell that we are studying varies 
less than a given percent, and during -at least- an interval of 
time equal to once the time of turnaround. In practice, what 
we do is to choose the moment immediately after the time 
of first shell-crossing, and we see if there is no a variation 
in radius larger than the maximum variation that we want 
to impose as our criterion. It must be in this way during, at 
least, a time of turnaround from this moment onward. The 
value of reference that we take to measure the variations in 
radius is the value of the radius at the initial moment of this 
interval. If the stated percentage of variation of the radius 
is exceeded in any time within this interval, then we move 



ahead in time until we find an interval of time where the 
criterion is satisfied. The first moment at which this occurs 
is our time of stabilization, and the value of the radius at the 
time of stabilization is taken as our radius of stabilization. 

This is illustrated in Figure [J] where the moments of 
stabilization are shown for a 5% and 10% of allowed vari- 
ation of radius, for a particular cosmology (Sim ~ 0.3, 
SIa ~ 0.7), a value of Mfrac, and for different virial masses. 
As one can see, in the case of 10% the stabilization is reached 
inmediately after the first shell-crossing, but the stabiliza- 
tion according to only 5% of allowed variation of the radius 
is not reached until roughly half of a time unit later (in 
units of the turnaround time). Moreover, this stabilization 
remains during more than once the time of turnaround from 
that moment onward (which is the minimum required by our 
criterion). Only at late times, where the evolution is dom- 
inated by numerical noise, the stabilization becomes worse 
than 5%. However, the most important conclusion is that the 
stabilization is not reached in any case for a radius that is 
1/2 the radius of turnaround, i.e. the value that corresponds 
to the collapse factor assumed in the standard derivation 
of (S;,t,ir and 5vir- This fact constitutes another proof that 
tell us how inappropriate are the current assumptions done 
about the virialization. 

A detailed study was done for different values of M„ir, 
Mfrac and different cosmologies. A summary of this study 
can be found in Appendix B, Tables IBll and IB2I In these 
tables we present the values for the density contrasts that 
we obtain for two different degrees of stabilization (5% and 
10%), varying M^ir, Mfrac and for the Einstein-deSitter and 
S7a 7^ cosmologies. 

We must note that the stabilization we have found is nu- 
merically robust, in the sense that it is independent of the 
integration parameters. In fact, it is in order to be certain 
that this result is not an artifact of approximation used for 
the dynamics, that we have used the direct numerical inte- 
gration of the equations of motion. As we shall show in future 
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Figure 4. Stabilization at 5% (circles) and 10% (crosses) for a 
particular cosmology (f2m = 0.3, Qa = 0.7) and Mj^^c = 0.5 
for three different virial masses. From down to top, the curves 
correspond to M„ir = 6.5 X 10'^'^h~^MQ M^ir = 3 X IO'^^/i'^Mq 
M^ir = 5 X l(!i^'^h~^ Mq. Again, the horizontal dashed-line cor- 
responds to half the turnaround radius, and the vertical one the 
time of collapse, i.e. twice the turnaround time. 



work, non-radial orbits and triaxiality may modify the form 
of the stabilization, but here we show that shell-crossing it- 
self leads to a stabilization, although in the stabilization of 
actual objects other dynamical factors are probably domi- 
nant. The origin of the stabilization lies on the fact that, 
as soon as a given shell is crossed, an outward particle flow 
develops which very nearly cancels the initial inward flow 
(at the given shell just before shell crossing) . It may also be 
understood as resulting from the radial "pressure" which ap- 
pears after shell crossing, that counterbalances gravity. Ex- 
plaining this fact in the simplest possible terms is a rather 
interesting initiative that we shall pursue. However, our in- 
tention in this work with respect to this point is to establish 
the fact that cold collapse (i.e. with no velocity dispersions) 
of a spherically symmetric cloud with a declining profile does 
stabilize. 



4.2 Virialization 

Concerning the virialization, what we did was to calculate 
the kinetic and potential energies related to the shell that 
we want to study. Then, we estimate the degree of agree- 
ment with respect to that given by the virial theorem (i.e. 
U -\- 2K = 0, where U and K are the potential and kinetic 
energies related to the shell under study) by means of the 
quantity: 



VIR- 



\U + 2K\ 
2K 



(26) 



It must be noted that, if the virial theorem was exactly 
satisfied by the shell at some moment of its evolution, this 
quantity should be at that time equal to zero. 

The algorithm to define the moment and radius of virial- 
ization is the same as already described for the stabilization. 



When we find an interval where the degree of virialization 
that we want to impose is satisfied in every moment inside 
this interval, then we define our time and radius of virial- 
ization as those ones corresponding to the beginning of the 
interval considered. Again, as in the case of the stabiliza- 
tion, a detailed study was done for different values of M^ir, 
Mfrac and different cosmologies. We summarize the results 
found in Tables lB3l and lB4l of Appendix B. In these tables we 
present the values for the density contrasts that we obtain 
for two different degrees of virialization (15% and 25%). 

There is one issue related to Tables IB3I and IB4I that 
it is worth mentioning. In those tables, a percent of 15% 
and 25% was set to look for virialization. This was done in 
this way because we noticed that below these percents it is 
impossible to reach virialization in most cases according to 
our criterion. A smaller percent means a degree of virializa- 
tion too strong to be satisfied. Nevertheless, it seems that, 
when the virialization is reached using these high percents, 
the halo has actually reached the virialization. This may 
be inferred from the fact that the individual percents that 
we measure in every moment within the interval considered 
does not decrease monotonically from the beginning of the 
interval to its end. In fact, what one obtains is a small fluctu- 
ation around (and near) the high percent that was imposed 
to find the virialization. That is, the degree of virialization 
does not vary substantially within the whole interval, only 
small fluctuations are found. There is a possible explanation 
to the fact that we have actually reached virialization but 
the degree of virialization that we find according to our def- 
inition is still above 10% or more. Until now, we have used 
the standard virial theorem, that only involves the kinetic 
and potential energies. However, because of the fact that 
we are treating with a non-isolated system, with shells of 
matter continuously going in and going out from the system 
that we are considering as our halo, we should include in the 
theorem another extra term. This term would be related to 
the pressures involved in the system, and surely may be the 
explanation and the cause of this "residual" percent that we 
obtain in all the cases. 



4.3 Comparison between stabilization and 
virialization: general considerations 

A summary of our results concerning the degree of both 
stabilization and virialization for a given shell can be found 
in Tables |4] and [5] In these tables, we show the degree of 
stabilization (STA) and virialization (VIR) reached for two 
known moments of evolution, A180 and A340 (see section 
3), for different values of M^ir, Mfrac and different cosmolo- 
gies. The degree of virialization, VIR, is calculated using 
Eq. (|26|l for moments A180 and A340. Concerning the degree 
of stabilization, it was estimated by calculating the following 
factor: 



STA = 



rviR — rsTA 



(27) 



rsTA 

where tvir is the value of the radius when 5 = 180 
(A180) or when 5 = 340 (A340), and rsTA is the value of 
the radius at the moment of stabilization, this one calculated 
according to the method described in section and using a 
value of 10% for the allowed variation in radius. This factor 
STA can be understood as a factor that measures the relative 
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Table 4. Values of VIR and STA, calculated according to equations H26|l and l|27| l. for two moments of evolution A180 and A340, and 
for different values of Mj^ac ^^id cosmologies. A virial mass of M„i^ = 3 X lQi^'^h~^ Mq was used in all the cases. (3 stands for the 
ratio between Q,m and Ha at the corresponding moment. See text for details. 



Einstein-deSitter {Vim = 1, f^A = 0) 





Mfrac 


= 0.2 


Mr 
frac 


= 0.5 


frac 


= 0.8 


Mfrac 


= 1.0 


frac 


= 1.3 


Moment 


VIR 


STA 


VIR 


STA 


VIR 


STA 


VIR 


STA 


VIR 


STA 


A180 


6.24 


0.57 


2.96 


0.43 


1.70 


0.33 


1.25 


0.25 


0.83 


0.19 


A340 


5.31 


0.47 


2.52 


0.31 


1.47 


0.19 


1.09 


0.09 


0.31 


0.09 



Qa^O 





Mfrac — 


0.2 


Mfrac — 


0.5 


Mfrac = 


0.8 


I^Ifrac = 


1.0 


frac — 


1.3 


Moment 


VIR STA 




VIR STA 


P 


VIR STA 


/9 


VIR STA 


/3 


VIR STA 


P 


A180 


4.23 0.49 


1.959 


2.69 0.35 


1.478 


1.70 0.23 


1.162 


1.43 0.19 


0.999 


1.05 0.13 


0.801 


A340 


3.88 0.39 


1.843 


2.41 0.21 


1.383 


1.47 0.07 


1.085 


1.29 0.05 


0.930 


0.95 0.07 


0.742 



contraction of the radius at times A180 or A340 with respect 
to that found at the moment of stabiHzation. 

Both tables provide useful results to extract important 
conclusions. In first place, one can see that the points A180 
and A340 are really far from the stabilization and also from 
the virialization, although these points are the preferred val- 
ues for the moment of virialization in most of the works 
found in the literature. In fact, in all the cases the first 
shell- crossing occurs even after A180 and A340, as pointed 
in section O so it is not possible that the shell has reached 
virial equilibrium or simply stabilization in radius at that 
moment. Also Lemson (1995) found similar results using N- 
body simulations, i.e. the equilibrium is reached in a longer 
time respect to that predicted by the standard SIM. In most 
of the cases, both stabilization and virialization were ob- 
tained far from the value 5i = 1.686 or 5i = 1.676, the 
preferred values of 5!,„ir for the Einstein-deSitter and the 
f2m = 0.3, r^A = 0.7 cosmology respectively (see Appendix 
B, Tables lBl1 to lB4|l . Concerning the associated values of the 
actual density contrast, 5, we find very high values in most 
cases. It should be noted, however, that it is expected that 
these values decrease substantially when we include angular 
momentum and velocity dispersion in the formalism. There- 
fore, the values showed in these tables are totally related to 
the isolated effect of shell- crossing. 

There are also other issues that could be interesting to 
stress, and that should be explored in more detail in a future 
work: 

(i) Concerning A180 and A340, the degree of both virial- 
ization and stabilization are better for larger values of M-uir- 

(ii) Concerning A180 and A340, the degree of both viri- 
alization and stabilization are also better for larger values 

of Mfrac- 

(iii) Concerning A180 and A340, the degree of both viri- 
alization and stabilization are worse in the Oa cosmol- 
ogy- 

(iv) The moment of stabilization is reached earlier in the 
JIa 7^ cosmology. 

(v) The moment of virialization is reached earlier for 
smaller values of M^ir- 

(vi) The moment of virialization is reached earlier for 



larger values of Mfrac- It is worth mentioning that Lemson 
(1995) found the same from his simulations, i.e. the inner 
shells reach equilibrium later. 

(vii) The moment of virialization is reached later in the 
f^A 7^ cosmology. 

(viii) It seems that there is no coincidence in time between 
virialization and stabilization, although it should substan- 
tially depend on the percent that we impose to define both 
concepts. 

(ix) For all values of Mfrac and Mvir, the stabilization in 
radius occurs at a fraction of the turnaround radius that is 
different from that given by the collapse factor (1/2 in the 
Einstein-deSitter case). The same is valid for the f^A 7^ 
cosmology. 



5 SUMMARY AND FUTURE WORK 

In this work we have studied the effect of shell-crossing in 
the formation and subsequent evolution of dark matter ha- 
los. To do that, we have used the spherical collapse model, 
which has been widely used in the literature for more than 
thirty years to manage and solve questions related to these 
processes. Despite of the large amount of works that have 
used this model or many others that have improved it by in- 
troducing in the formalism more and more complex consid- 
erations, only a few of them have included the effect of shell- 
crossing. Moreover, most of these works have managed this 
effect analytically using the adiabatic invariant as a good 
approximation . 

Here we handle the effect of shell-crossing numerically. 
This allows us to study individually any shell of matter in- 
volved in the process of formation of the halo. Doing so, 
we can extract multiple conclusions about the way in which 
this process occurs, like the relation between the linear and 
actual density contrasts, the process of stabilization of a 
shell of matter, the virialization, etc. Most of these issues 
have been treated in the present work to a greater or lesser 
extent, although the main goal have been always the devel- 
oping of an adequate framework - named as Spherical Shell 
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Table 5. Values of VIR and STA, calculated according to equations H26| l and l|27| l. for two moments of evolution A180 and A340, and 
for three different virial masses and two cosmologies. A value of Mf^ac = 0.5 was used in all the cases. /3 stands for the ratio between 
Qrn and £7a at the corresponding moment. See text for details. 



Einstein-deSitter {Qm = 1, f^A = 0) 







= 6.5 X lOio?t-iA/0 




= 3 X W^h-^MQ 


A4 


= 5 X lO^'^h-'^MQ 


Moment 


VIR 


STA 


VIR 


STA 


VIR 


STA 


A180 


2.75 


0.51 


2.96 


0.43 


2.00 


0.23 


A340 


2.16 


0.41 


2.52 


0.31 


1.92 


0.05 



Qa^O 







= 6.5 X 10 






= 3 X 10^2 


h-^MQ 




= 5 X 10" 


h-^MQ 


Moment 


VIR 


STA 




VIR 


STA 


/3 


VIR 


STA 




A180 


3.45 


0.43 


1.340 


2.69 


0.35 


1.478 


1.16 


0.17 


1.849 


A340 


2.80 


0.31 


1.254 


2.41 


0.21 


1.383 


1.12 


0.05 


1.737 



Tracker (SST) - in which we can study in depth the sheU 
crossing and also other secondary effects. 

It is possible to summarize the main conclusions of this 
work as follows: 

(i) The SST framework is adequate to tackle the effect of 
shell-crossing in a way that allow us to extract exact results 
for different issues related to the evolution of the halo: the 
way that the radius of a given shell evolves with time, the 
relation between the linear and actual density contrasts, the 
stabilization, the virialization, etc. 

(ii) The shell-crossing by itself is able to produce stabi- 
lization and virialization. Nevertheless, for the moment, it is 
not possible to obtain the exact values of the linear and ac- 
tual density contrasts related to both moments of evolution. 
It is necessary to take into account also other important 
effects, such as angular momentum, velocity dispersions or 
triaxiality. 

(iii) Concerning the relation between the linear and ac- 
tual density contrasts, the dependence of this relation with 
the cosmology is very small and practically negligible. This 
conclusion is contrary to most of previous works, which find 
in general a large dependence with cosmology. However, the 
dependence with the virial mass or the fraction of virial mass 
that we consider, is large. 

(iv) Neither stabilization nor virialization are reached in 
a time according to that given by the common assumptions 
related to the collapse factor and the time of virialization. In 
all the cases, we find that both stabilization and virialization 
occur at later times. 

(v) The values typically used in the literature for 5i,vir 
and Svir seem to be clearly inadequate and incorrect, and 
are based on not very solid assumptions. In this work, new 
values of Si^vir and S^ir are presented, but only taking into 
account the effect of shell-crossing. It will be necessary to 
include in our framework other effects also relevant to be 
able to provide useful and final values for 5i,vir and Spir- 
it is worth to emphasize that this work constitutes only 

a first step in our attempt to obtain exact and precise pre- 
dictions related to the formation and evolution of dark mat- 
ter halos. In a future work we plan to include in the SST 



framework other important effects that it will be absolutely 
necessary. In particular, including the angular momentum 
and velocity dispersion will be the next step. Furthermore, 
in parallel, we will implement a more sophisticated initial 
density profile than that used in this work, which fits better 
that found in the simulations and could change the results 
presented here only slightly. It is also in our mind to use 
cosmological N-body simulations, since comparison between 
both analytical and simulation studies will be, by sure, cru- 
cial to reach a better and deeper understanding of the pro- 
cesses involved in the formation and evolution of dark matter 
halos. 
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APPENDIX A: THE FORMALISM IN THE 
KCDM COSMOLOGY 

If we are in a A 7^ cosmology, we need to introduce some 
changes for the initial conditions and in the expressions for 
the evolution of the spherical perturbation, although the for- 
malism and the algorithm are essentially the same as pre- 
sented in section 2, i.e. the SST framework. 

The equation for the initial radii of the shells is the 
same as given by Ea. (|18|) . but for the velocities the correct 
expression, instead of Eq. (|19|l . is now: 



1 



ffli b{a) 



3 1+5(5?(g(j))) D{a) 



dSi(S) I 



X n(j) 3 + In 



where 



(Al) 
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dD{a) 
da 



and, 



a{t) 



PS 



(ft ^ + vi+^)* - (ft ^ + \n + W^y 



1 



1 + 2. 

D{a) = 



1 



/(«) 



2 a f{a] 

1 +fim 



/ (a) da 



- Q,A{a 



(A2) 



Note that /3i is simply the parameter /3o = f2m/f^A but 
referred to the initial time. 

Concerning the initial density profile, now the initial 
Unear density contrast is essentially the same as given by 
Ea. (|15p but now the rescaling factor, yqri") is here replaced 



by D{a=\) ' which gives: 



SiiqU)) 



(A3) 



D{a = 1) 

where Si{q{j)) is the linear profile given by Eq. (|12p . 

For the evolution, equations (I21|l and (I22p are still valid 
for the radius and the velocity, but to compute the linear 
and actual density contrast, now we have to include: 



[1 + S{5i{q{m 



ai 



(A4) 



(A5) 



where D{a) and a{t) axe the growing and scale factor 
respectively, as defined in Ea. HA2p . and ai denotes the scale 
factor at initial time, given also in Eq. (|A2|l . 

To recompute the enclosed mass at each time step, it 
is also necessary to take into account the new cosmology, 
once we have calculated M{j, t) in first place according to 
Eq.((25l), i.e.: 



(A6) 



M(j,t)(A / 0) = M[j,t)[K = 0) - -r{j,tf 

Pi 



APPENDIX B: RESULTS OBTAINED FOR 
STABILIZATION AND VIRIALIZATION 

In Tables IBll to IB4I the linear and actual density contrasts 
that we obtain concerning the moments of stabilization and 
virialization are shown. Tables IbTI and [B2l refer to the sta- 
bilization whereas Tables [B3l and lB4l are related to the viri- 
alization. In both cases, the moments of stabilization (STA) 
and virialization (VIR) were matched following the criteria 
given in Section|4l In these Tables, linear and actual density 
contrast are shown for different values of virial mass, Mvir, 
fraction of virial mass, Mfrac, and for Einstein-deSitter and 
Q.A 7^ cosmologies. 
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Table Bl. Linear and actual density contrasts for two different values (0.05 and 0.10) of STA, this one defined according to expression 
II27I I. for different values of Mf^ac and two cosmologies. A virial mass of M^i,- = 3 X 10^'^h~^ Mq was used in all the cases. /3 stands for 
the ratio between Qrn and Ha at the corresponding moment. See Section |4] for details. 



Einstein-deSitter {Qm = 1, f^A = 0) 







Mfrac 


= 0.2 


Mfrac = 0.5 


Mfrac = 0-8 


Mfrac = 1.0 


Mfrac 


= 1.3 






STA 


Si 


S 


Si s 


Si s 


Si s 


Si 


5 






0.05 




















0.10 


1.98 


3993 


2.00 1789 


2.05 1223 


1.89 700 


1.90 


537 














Qa^O 












Mfrac = 0.2 


Mf 


7'ac — 0-5 


M frac = OS 


Mfrac = 


1.0 


frac — 


1.3 


STA 


Si s 


/9 


Si 


s p 


Si S 13 


Si s 


/3 


Si s 


/3 


0.05 


1.91 3306 


0.959 


1.85 


1298 0.813 


1.86 897 0.577 


1.86 741 


0.483 


1.86 621 


0.362 


0.10 


1.67 1496 


1.697 


1.67 


703 1.269 


1.67 463 0.977 


1.67 390 


0.829 


1.67 324 


0.663 



Table B2. Linear and actual density contrasts for two different values (0.05 and 0.10) of STA, this one defined according to expression 
II27I I. for three different virial masses and two cosmologies. A value of Mfrac = 0.5 was used in all the cases. (3 stands for the ratio 
between and Ha at the corresponding moment. See Section |4] for details. 



Einstein-deSitter (H^ = 1, S^A = 0) 



STA 


Si 


= 6.5 X lO^^ft-iM© 
6 


M^ir = 3 X lO^2/i"^A/0 
Si s 


Si 


= 5 X lO"/i-^M0 
5 


0.05 
0.10 






2.00 1789 


1.87 


610 


Ha ^0 


STA 


s, 


= 6.5 X I0^°h-^MQ 
S f3 


M^ir = 3 X IO^^/i^^Mq 
Si S (3 


Si 


= 5 X lO"/i-Uf0 
S (3 


0.05 
0.10 


1.93 
1.67 


2727 0.569 
1105 1.137 


1.85 1298 0.813 
1.67 703 1.269 


1.89 
1.67 


680 0.944 
352 1.602 



Table B3. Linear and actual density contrasts for two different values (0.15 and 0.25) of VIR, this one defined according to expression 
H26| l. for different values of Mfrac a^nd two cosmologies. A virial mass of M^ir = 3 X lO^^h-iM© was used in all the cases. (3 stands for 
the ratio between f2m and Ha at the corresponding moment. See Section |4] for details. 



Einstein-deSitter (^m = 1, f^A = 0) 
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Si S 13 


Mf 

Si 
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S 13 
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Si s 


.8 

/9 


Mfrac = 

Si s 


= 1.0 

/3 




Mfrac = 1.3 
Si S f3 


0.15 
0.25 


2.65 23382 0.119 


2.47 


8359 0.119 


2.33 4939 


0.119 


2.25 3119 


0.119 




2.14 2027 0.119 
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Table B4. Linear and actual density contrasts for two different values (0.15 and 0.25) of VIR, this one defined according to expression 
1261 1. for three different virial masses and two cosmologies. A value of Mf^ac = 0-5 wS'S used in all the cases. /3 stands for the ratio 
between Qm and Q\ at the corresponding moment. See text for details. 



Einstein-deSitter {Qm = 1, f^A = 0) 
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Si 


= 6.5 X lO'^^h^'^MQ 
5 


M^ir = 3 X IO^^/i-iMq 
5i 5 


Si 


= 5 X lOi4/i-iAf0 
5 


0.15 
0.25 


1.92 




2132 


2.20 2513 
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2.91 
2.72 




3886 
2378 




VIR 


Si 


= 6.5 X lO"^/i"^M0 

s p 


M^ir = 3 X W^^h-'^MQ 
5i S /3 


A'''Iy i f 

Si 


= 5 X lO"/i"^Af0 
5 /3 


0.15 
0.25 


2.19 
2.08 


6056 
4141 


0.268 
0.376 


2.47 8359 0.119 


2.53 
2.42 


3931 
3177 


0.161 
0.237 
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